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ABSTRACT 


I We  establish  conditions  for  the  appearance  of 
nonlinear  -^cusped  waves-^in  density  layered  shear  flows 
by  means  of  computer  simulations  with  a model  set  of' 
initial  mean  density  and  shear  ve]ocjty  profiles.  These 
waves  occur  when  a region  with  gradient  Richardson 
number  R*  > 1/4  is  initially  confined  between  regions  where 

R*  < 1/4.  In  other  cases  a closed-circulation,^ -^cat ' s eye1®  1 

pattern  develops.  Buoyancy  forces  dominate  in  the  nonlinear 
evolution  of  cusped  waves,  whereas  Reynolds  shear  stresses 
dominate  in  the  development  of  cat's  eyes.  Cusped  waves 
continue  to  extract  energy  from  the  mean  shear  approxi- 
mately linearly  in  time,  whereas  cat's  eyes  have  bounded 
energy  content.  Indirect  evidence  is  presented  for  the 
existence  of  cusped  wave  Instabilities  in  the  ocean. 


■psttg t*0 


1 . Int.roduct,  inn 

Extensive  horizontal  layers  are  a persistent 
feature  of  stably  stratified  flows.  These  have  been  observed 
in  vertical  profiles  of  temperature  and  salinity  in  the  oceans 
(e.g.,  Elliott,  et_  aJ . , 197*0  and  in  the  temperature  structure 
of  the  atmospheric  boundary  layer  during  night-time  radiation 
inversion  (Turner,  1973)-  Photographs  by  divers  in  the 
seasonal  thermocline  of  the  Mediterranean  (Woods,  1968)  have 
shown  that  a long  wavelength  internal  wave  propagating  on  a 
thin  interfacial  "sheet"  between  different  constant  temperature 
and  density  layers  can  lead  to  shear  instabilities  of  the 
Kelvin-Helmholtz  type.  Examples  of  shear-induced  insta- 
bilities have  also  been  observed  in  the  atmosphere  with 
radar  sounders  (Gossard  and  Hooke,  1975).  Similar  examples 
of  this  instability  have  been  studied  in  the  laboratory  by 
Thorpe  (1973),  using  a density  layered  fluid  tube. 

The  nonlinear  development  of  shear  - induced 
instabilities  in  the  oceans  provides  a mechanism  for  enhanced 
vertical  transport  of  heat  and  momentum.  Shear  instabilities 
may  therefore  be  an  important  component  in  thermocline 
formation  (Munk,  1972).  Woods  and  Wiley  (1972)  have  suggested 
that  many  of  the  observed  details  of  thermocline  temperature 
structure  can  be  attributed  to  shear  instabilities  caused 
by  long  internal  waves  propagating  on  sheets. 


According  to  these  authors,  the  turbulent  billows  which  form 
entrain  fluid  from  either  side  of  the  interfacial  region 
until  the  sheet  thickness  grows  to  several  times  its  original 
value.  Internal  wave  shear  eventually  causes  neighboring 
billows  to  overlap,  providing  a possible  mechanism  for 
replacement  of  the  original  interface  by  two  sheets.  Subse- 
quent internal  waves  lead  to  further  splitting,  resulting  in 
vertical  transport.  Garrett  and  Munk  (1972)  have  used  these 
ideas  to  estimate  the  depth  dependence  of  vertical  and 
horizontal  diffusivities  in  the  ocean. 

The  experiments  of  Thorpe  (1973)  indicate  that 
the  sheet-splitting  mechanism  of  Woods  and  Wiley  requires 
further  investigation.  The  vertically  spreading  turbulent 
regions  observed  by  Thorpe  show  no  evidence  for  large  gradients 
In  density  or  horizontal  velocity  at  their  edges.  These 
experiments  also  disagree  with  a qualitative  description  of 
sheet  spreading  advanced  by  Businger  (1969).  The  latter 
argues  that  spreading  ceases  when  the  mean  gradient  Richardson 
number  of  the  turbulent  region  rises  to  unity.  Thorpe, 
however,  measures  fully  developed  mean  gradient  Richardson 
numbers  closer  to  0.3* 

The  present  paper  discusses  the  results  of 
extensive  two-dimensional  numerical  simulations  which  address 
the  above  issues.  The  fully  developed  states  which  result  from 
shear  instabilities  on  thermal  sheets  in  the  ocean  possess 
several  features  which  make  them  attractive  for  numerical 
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simulation.  They  are  first  of  all  characterized  by  moderate 
Reynolds  numbers  Rg  (based  on  sheet  thickness  and  mean  shear 
velocity)  of  order  1000  (Thorpe,  1973),  corresponding  to 
microscale  Reynolds  numbers  of  order  100.  In  contrast  to 
most  other  interesting  turbulent  geophysical  flows,  numerical 
simulations  with  Reynolds  numbers  of  this  order  are  possible 
without  introducing  phenomenological  treatments  of  Reynolds 
stress  terms  (Fox  and  Lilly,  1972;  Case,  et  al . . 1973). 
Buoyancy  forces  also  act  to  suppress  vertical  turbulent 
motions  and  to  couple  energy  into  large  -scale  internal  waves. 
Finally,  the  initial  shear  instability  and  many  features  of 
the  fully  developed  billows  are  two-dimensional  effects. 

Calculations  are  presented  for  Prandtl  numbers 
Pr  = 10,  relevant  to  thermal  layers  in  the  ocean,  and  for 
Reynolds  numbers  Re  <_  200,  slightly  lower  than  typical 
billow  mixing  events  In  the  ocean  and  in  Thorpe's  experi- 
ments. However,  results  obtained  at  these  Reynolds 
numbers  are  sufficiently  similar  to  observed  phenomena  to 
be  of  considerable  interest.  In  particular,  the  origin  of 
"cusped  waves"  is  explained.  These  waves  have  been  observed 
in  a variety  of  turbulent  mixing  experiments  (Ellison  and 
Turner,  1959;  Lofquist,  I960;  Thorpe,  1968;  Moore  and  Long, 
1971)  and  in  the  atmosphere  (Gossard  and  Hooke,  1975). 

Thorpe  (1973)  states  that  these  waves  appear  at  about  the 
time  of  the  onset  of  turbulence.  They  first  appear  (or  break) 
in  the  upper  layer  of  fluid,  but  at  later  times  they  are 
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observed  moving  In  opposite  directions  on  either  side  of 
the  interface. 

The  results  of  the  present  paper  confirm  a conjec- 
ture by  Thorpe  (1973)  that  cusped  waves  are  the  finite  ampli- 
tude result  of  shear  instabilities  associated  with  small 
gradient  Richardson  numbers  at  the  sheet  edges  and  large 
Richardson  numbers  >_  1/4  within  the  sheet.  Calculations  show 
that  this  instability  is  asymmetric  and  first  appears  in  the 
upper  portion  of  the  flow.  In  other  cases,  where  the 
Richardson  number  is  less  than  1/4  within  the  sheet,  the  insta- 
bility is  symmetric  and  develops  into  a familiar  "cat's  eye" 
circulation  pattern.  The  latter  situation  has  been  studied 
numerically  by  Patnaik,  et^  al.  (1976),  for  Pr  < 1 . i 

Our  calculations  also  reveal  that  buoyancy  forces 
dominate  Reynolds  stresses  In  the  evolution  of  cusped  waves, 
so  that  these  are  a truly  buoyant  phenomenon.  The  reverse 
is  true  in  the  development  of  the  cat's  eye  circulation  pattern, 
which  is  a basically  nonbuoyant  phenomenon.  Moreover,  cusped 
waves  continue  to  extract  energy  from  the  mean  shear  and  to 
increase  the  sheet  thickness  at  very  long  times;  the  cat's  eye 
pattern  grows  to  a maximum  vertical  extent  and  energy  content 
and  then  slowly  decays  by  viscosity. 

The  possible  importance  of  the  cusped  wave 
instability  in  the  ocean  is  examined  in  a concluding  section. 
Using  a model  thermocline,  similar  to  that  reported  by 
Woods  (1968)  in  the  summer  thermocline  near  Malta,  It  Is 


L 


shown  that  the  velocity  and  thermal  structure  of  an  internal 
wave  modified  by  the  presence  of  a thermal  sheet  is,  in  many 
cases,  more  favorable  to  the  formation  of  cusped  waves  than 
cat's  eyes. 


2 . Linear  Stability 

We  first  present  relevant  properties  of  shear 
induced  instabilities  due  to  a thin  horizontal  sheet  of 
enhanced  thermal  gradient.  We  consider  unperturbed  horizontal 
velocity  and  density  profiles 


u(z)  = UQerf  ^ i /?  z/6^ , 

p(z)  = poexp^-0/fi  erf^j-  S*  Hz/ 6^  , 


(2.1) 


(2.2) 


where  z (-<*>  < z < °°)  is  the  vertical  coordinate,  Uq/6  and 
Pq  are  the  horizontal  shear  and  density  at  the  sheet  center 
z = 0,  R is  the  ratio  of  half-widths  of  the  velocity  and 
density  profiles,  and  8 = R sinh”^^  Ap/Pq^.  Here  Ap,  the 
total  density  difference  across  the  sheet,  is  related  to  the 
thermal  difference  AT.  For  example,  in  water  near  20°C,  small 
density  and  thermal  variations  are  related  by 


--  = -0.057 


(2.3) 


where  the  mean  temperature  is  expressed  in  degrees  Kelvin. 
Hazel  (1972)  ha::  studied  similar  (hyperbolic  tangent)  profiles 
for  various  thickness  ratios  K. 
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In  the  following,  v/o  make  quantities  nondimensional 
with  appropriate  combinations  of  6,  Uq  and  Pq.  The  local  gradient 


Richardson  number  corresponding  to  (2.1)  and  (2.2)  becomes 


Ri(z)  = J exp  - jj-  (R^-2)z2 


(2.4) 


where 


J = Bg6/U‘ 


is  the  overall  Richardson  number  based  on  the  sheet  thickness 


6.  For  R < /2  , the  gradient  Richardson  number  has  a minimum 


at  the  sheet  center  and  profiles  (2.1)  and  (2.2)  are  stable 


for  J >_  -jj-  . Figure  1 shows  neutral  stability  curves  which 


we  have  computed  from  the  inviscid  Taylor  - Goldstein 
equation  in  the  Boussinesq  approximation 


(u-c)  = 0. 

\ dz^  / V dz  u - c / 


(2.5) 


Here  <J>  ( z)  is  the  vertical  profile  of  the  perturbation  stream 
function 


ip'(x,z,t)  = Re|<J>  ( z ) exp[ia(x-ct ) ]|  (2.6) 


and  a and  c = c^  + ici  are  a dimensionless  wave  number 


and  complex  phase  speed.  N(z)  is  a dimensionless  mean 
Brunt-Vaisala  frequency 


n2<z>  ' * ^ Si  ln  » - 


(2.7) 


whore  Fp  = i.;  u Froudc  number. 


The  growth  rates  ae^  of  the  most  unstable  disturbances  are  shown 

(dashed  line)  for  R = 1 in  the  upper  half  of  Figure  1.  The  lower 
part  of  the  figure  shows  neutral  stability  curves  (where  c = 0) 
as  a function  of  R.  For  a < /2  the  neutral  stability  curves 
bound  the  region  of  instability  in  the  a - J plane. 

For  R > /2  (shear  profile  wider  than  density  profile) 
the  gradient  Richardson  number  approaches  zero  for  large  \z\ . 
Figure  2 shows  stability  boundaries  for  R = 3-5.  The  neutral 
stability  curve  bounds  the  region  of  stability  for  one  mode, 
but  another  mode  is  always  unstable  (Hazel,  1972).  Inside 
of  the  solid  curve  the  instability  is  a standing  wave 
(c  = 0)  and  outside  it  is  progressive.  The  figure  also  shows 
that  the  locus  of  the  most  unstable  mode  moves  to  shorter 
wavelengths  as  J is  Increased  above  1/4. 

Figure  3 shows  real  and  imaginary  parts  of  the 
complex  stream  eigenfunction  <p(z)  normalized  so  that  <J>(0)  = 1. 
The  most  unstable  <p(z)  is  shown  for  various  R and  J.  Case  (b) 
is  typical  for  R < /2  and  does  not  differ  significantly  from 
the  J=0  unstratified  case  (a).  For  case  (c),  the  entire 
flow  satisfies  R1  c jj-  and  the  eigenfunctions  are  relatively 
broad.  Case  (d)  is  a propagating  mode  (c  = 0.330)  with 
R = 3-5  and  J = 0.4.  The  most  unstable  perturbation  is 
asymmetric  about  z = 0.  As  we  shall  see,  this  asymmetry 
leads  to  nonlinear  cusped  waves  in  the  upper  half  of  the 
shear  region. 
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3 . Numerical  He  horn* 

3 • 1 Bousslnesq  Equations 
We  assume  that  local  density  variations  are  small 
with  respect  to  pQ  and  are  linearly  related  to  thermal 
variations,  as  in  (2.3).  We  shall  consider  the  nonlinear, 
two-dimensional,  late  time  evolution  of  small  perturbations 
of  profiles  (2.1)  and  (2.2)  according  to  the  Boussinesq 
equations 


9tu  + V-(uou)  = - Vp  - F~1p£  + R”1V^u , (3.1) 

3tP  + V • ( pu ) = Pg^p,  (3.2) 

V-u  = 0.  (3.3) 

Here  o denotes  a dyadic  tensor  product, 

V = Ox,3z),  (3.4) 

u(x,t)  = [u(x,z,t) ,w(x,z,t)]  (3.5) 


is  the  two-dimensional  velocity  in  a vertical  x - z plane, 
p is  the  nondimensional  pressure,  z is  the  unit  normal  in 
the  vertical  direction,  Rg  is  the  Reynolds  number 

U 6 

Re  " -T-  * (3-6) 

v is  the  kinematic  viscosity,  ?e  is  the  Peclet  number 

re  = RoPi-  * (3.7) 


P is  the  Praridtl  number 

r 


Pr  = v/k 


and  k is  the  thermal  diffusion  coefficient.  In  the  following 

we  use  values  = 10  and  F = 0.001,  representative  of  the 
r r 

ocean  (Woods,  1968). 

3 . 2 Boundary  Conditions 

We  use  periodic  boundary  conditions  for  0 <_  x <_  2ir/a 
in  the  horizontal  direction.  We  confine  the  flow  between 
horizontal  planes  - ^ £ z <_  ^ in  the  vertical  direction  and 
impose  finite  boundary  conditions 


3zu  = 0,  w = 0,  p = p+  or  p_  (z  = ± D), 


(3.9) 


where  p+  and  p_  are  constant  density  values  at  the  upper  and 
lower  boundaries,  as  determined  by  (2.2).  Conditions  (3-9) 
guarantee  that  vorticity 


u = -3  u + 5 w 

Z X 


(3.10) 


is  not  produced  at  the  upper  and  lower  boundaries.  The 
pressure  boundary  condition  corresponding  to  (3*9)  follows 
from  (3.1)  as 

3zp  = -F~^p  = constant  |z  = * ^ d|  . (3-11) 

3 . 3 Spatial  Discretization 
We  use  a numerical  scheme  similar  to  that  described 


by  Zubusky  and  D<  <-m  (1971).  Accordingly , discretized 


u 


staggered  veloc  1 tier,  and  dour.  1 ties  u , v;  and 

m + 2,n  m,m  + ^ 

P 2 are  defined  on  lattice  cell  sides  and  pressures 
m,n  + p 


m 


>n  are  defined  at  cell  centers  for  m = 1 , . . . ,M  and 


n - 1 , . . . ,N.  Equations  (3.1)  and  (3-3)  are  replaced  by 


(3.12; 


6 

x 


u 

m 


1 

2 


0 w 

z 1 

m,n  - - 


0, 


(3.13) 


with  similar  equations  for  u l and  p . Here 

^ 2^  . n m,n  + 

anc^  are  one~sided,  forward  difference  operators  in  the 

x and  z directions,  respectively . Quantities  such  as  w a; 

m , n 

appropriately  centered  linear  averages 


w 

m,n 


w 

m, 


(3-14) 


etc.  In  the  absence  of  density  variations,  the  above  spatial 
discretization  guarantees  conservation  of  momentum  and 
energy  (for  v = 0) . 
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3 • ^ Temporal  Dlscr^t  1 sat  ion 

We  use  a split  -step  approach  to  advance  in  time. 

Auxiliary  velocities  u ^ and  w n are  first  computed 

m + 2>n  m,n  + |- 

using  a time  - centered , three  time  - level  " leapf rog"  scheme,  which 
omits  the  pressure  and  viscosity  terms.  Time  - centered  pres- 
sures are  then  computed  from  the  Poisson  equation 


* 

6 v 

z 

m ,n 


(3.15) 


with  horizontal  boundary  condition  (3.11),  where  At  is  the 
time  step.  Equation  (3-15)  is  solved  exactly  and  noniteratively 
using  cyclic  reduction  methods  (Buneman,  1969;  Schumann  and 
Sweet,  1976).  New  auxiliary  velocities  are  then  computed 
from 


**  * 

u 1 = u 1 
m + ^jn  m + 2-,n 


2At6  p 

x m ,n. 


(3.16) 
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Fur  largo  HuynoJda  numbers,  the  above  scheme  Is 

2 2 2 

second  -order  accurate  0(Ax  +Az  +At  ).  In  our  calculations, 
this  scheme  remains  stable  for  At  <_  min( Ax/u ,Az/w) , where 
Ax  and  Az  are  the  constant  x and  z lattice  spacings.  Moreover, 
it  retains  the  inviscid  time  reversibility  of  the  Navier- 
Stokes  equations  and  guarantees  incompressibility  (3.13)  to 
within  roundoff  error  of  the  computer. 

For  lower  accuracy  computers  It  is  preferable 
to  eliminate  the  large  static  pressure  by  replacing  p and  p 
in  (3.1)  with  the  perturbed  density  and  pressure 


P 


P 


P+  + P 
~ 


(3.18) 


P 


P 


P+  + P 
— 2 


(3-19) 


The  perturbed  pressure  then  satisfies  (3.15)  with  Neumann 
boundary  conditions 


(3-20) 


3 . 5 Initial  Conditions 

We  consider  only  the  growth  of  a single  perturbing 
mode,  using  initial  mean  profiles  (2.1),  (2.2)  and  the  complex 
eigenfunctions  p of  Figure  3 (altered  slightly  to  make  them 
vanish  at  the  boundaries  z = + D)  . From  (2.6)  and  the  per- 
turbed density  equation,  we  have 
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u(  x , s , 0 ) = u ( z ) + n3zRe| 'M-’)oxp(1ax)|  , (3-?l) 

(3.22) 


v(x , z , 0)  = - n3xRe  |<J>(z)exp(lax)| 


p(x,z,0)  = p(z)  + nR, 


<J>(z)  3 o(z) 


(z)  - c 


exp ( lax ) 


}• 


(3.23) 


Here  c is  the  complex  linear  phase  speed  and  n is  a small 
parameter  (0.1  or  less  in  our  computations). 

We  have  not  studied  nonlinear  interaction  between 
modes  (Zabusky  and  Deem,  1971)  or  subharmonic  generation 
(Patnaik,  et  al . , 1976)  in  the  present  paper,  although  these 
subjects  are  of  considerable  importance  in  observed  two 
dimensional  shear  flows  (Winant  and  Browand , 1973). 

In  succeeding  sections,  we  consider  separately 
cases  where  R < /2  and  R > /2,  where  R is  the  relative 
velocity  shear  thickness  in  (2.2).  Table  1 summarizes  the 
physical  and  numerical  parameter  ranges  which  we  have 
studied . 

4 . Nonlinear  Dynamics  for  R = 1 

Figures  4 and  5 show  the  nonlinear  development 
of  the  fastest  growing  mode  a = 0.43  for  R = 1 and  J = 0.1. 
Figure  4 gives  constant  density  contours  and  Figure  5 gives 
constant  vorticlty  contours  (solid)  and  streamlines  (dotted). 
The  flow  develops  into  a cat's  eye  pattern,  characterized 
by  closed  streamliner,  at  the  center  of  the  flow,  where  a 
significant  portion  of  fluid  becomes  confined.  After  t = 40 


! 
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the  pattern  la  fully  developed  and  begins  to  decay  due  to 
viscosity  and  thermal  diffusivity.  The  results  are  not 
significantly  altered  if  we  increase  D = 10  (Run  1)  to 
D = 16  (Run  2). 

Figure  6 shows  vertical  profiles  for  the 
horizontally  averaged  velocity  u(z),  density 
a( z)  = 1000* [p ( z ) -1 ] (Phillips,  1969)  and  gradient  Richardson 
number  R^z)  at  a late  time  t = 150  (Run  2).  Corresponding 
curves  at  t = 0 are  indicated  with  dashed  lines.  We  note 
that  R^(z)  approaches  zero  in  the  outer  flow  near  |z|  =5 
and  has  an  average  value  of  approximately  0.29  in  the  inter- 
vening flow  region. 

Figure  7 shows  the  time  dependent  behavior  of  the 
momentum  and  density  vertical  thicknesses,  normalized  with 
their  values  at  t = 0.  Here  Au  and  Ap  are  defined, 
respectively,  as  the  vertical  distances  between  z values  where 
u and  p come  to  within  90^  of  their  values  at  the  upper  and 
lower  edges  of  the  flow.  In  the  present  case,  where  R < /2, 
the  density  profile  spreads  by  less  than  a factor  of  two, 
whereas  the  momentum  profile  continues  to  grow  slowly, 
presumably  as  the  result  of  viscosity.  The  oscillations 
evident  in  the  figure  are  the  result  of  a time  periodic 
nutation  of  the  central  elliptical  streamline  pattern  in 
Figure  5 (Zabusky  and  Deem,  1971;  Patnaik,  e_t  aJ.  , 1976). 

Figure  8 shows  time  evolution  of  the  fluctuation 
kinetic  energy 


ifilVwa*- 


cr  = \ <n,?+w2>  , 


( 'i . i ) 


the  Reynolds  shear  stress 


S = - <u'w'*9zu'>  , 


(4.2) 


the  buoyancy  flux 


W = - F”1  < p ' w’  > , 


(4.3) 


and  the  dissipation  rate 


«=  R"1  <(3xu’)2  + Ozu’)2  + Oxw')2  + (9zw')2>  (4.4) 

for  Run  2.  Here,  we  divide  the  flow  variables  into  their  down- 
stream averaged  mean  and  fluctuating  parts  (e.g.,  u = u+u') 
and  we  use  braces  to  denote  averages  over  the  entire  flow. 

The  above  quantities  satisfy  the  fluctuation  kinetic  energy 
equation  (Phillips,  1969) 


^ Ef  = S + N - € . 


(4.5) 


so  that  S,  N and  e are  a measure  of  the  overall  contributions 
of  shear,  buoyancy  and  viscous  dissipation  in  the  balance 
of  kinetic  energy.  In  all  of  our  calculations  the  Reynolds 
number  Is  sufficiently  large  so  that  e is  of  secondary 
importance  (Figure  8). 

Figure  8 shows  a characteristic  feature  of  the 
cat's  eye  pattern,  namely  that  shear  stress  is  generally 
larger  than  buoyancy  flux. 
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5.  Nonlinear  DynamJcr.  for  R = 3-5 

For  R > /2  the  initial  gradient  Richardson  number 
falls  to  zero  for  large  |z|  and  has  a maximum  J at  z = 0. 

We  consider  cases  where  J < 1/4  and  J > 1/4  in  this  section. 

Figures  9 and  10  show  density,  vorticity  and 
stream  function  contours  for  the  nonlinear  development  of 
the  fastest  growing  mode  a = 0.38  with  R = 3-5  and  J = 0.1 
(Run  3).  Figure  11  shows  the  time  behavior  of  E^,,  S,  N 
and  e (Run  4).  As  one  would  expect,  the  instability  is  more 
energetic  than  for  R < /2  (Figure  8),  but  the  final  state 
is  similar  to  the  shear  stress  dominated  cat's  eye  pattern 
of  the  previous  section. 

Figures  12  and  13  show  that  cusped  waves  result 
for  R * 3*5  and  J = 0.4.  The  Reynolds  number  is  Rg  = 100 
(Run  7).  The  instability  is  nonsymmetric  and  propagating 


and  develops  preferentially  into  the  lighter  fluid  above 


the  central  shear  layer. 

Figure  14  shows  the  same  flow  at  very  long  times. 
The  density  contours  reveal  a secondary  wave  (arrows)  which 
moves  to  the  left  in  a direction  opposite  to  and  with  a speed 
equal  to  that  of  the  primary  cusped  wave.  Examination  of 
the  central  stream  function  contours  at  t = 160  and  t = 200 


shows  that  when  the  two  waves  are  in  close  proximity,  a closed 
circulation  region  forms. 

Figure  15  shows  that  at  higher  Reynolds  numbers 


200,  Run  9)  the  cusped  wave  break 


The  flow  is  again  characterised  by  a pair  of  nonlinear  waves 
moving  in  opposite  directions.  At  late  times  the  outer 
portions  of  the  lower,  secondary  wave  also  show  evidence  of 
breaking . 

/ 

We  have  performed  calculations  at  higher  Reynolds 
numbers  (Table  1),  but  the  present  numerical  resolution  is 
not  sufficient  to  resolve  all  the  fine  scale  features  which 
develop.  However,  we  do  call  attention  to  a trend  whereby 
the  upper  and  lower  vortex  regions  become  symmetric  at  late 
times . 

Figure  16  shows  the  time  behavior  of  E^.,  S,  N and  e 
for  Runs  7 and  9.  We  see  that  the  formation  of  cusped  waves  is 
a buoyancy  dominated  phenomenon  very  different  from  the  cat's 
eye  case.  Moreover,  cusped  waves  continue  to  extract  energy 
from  the  mean  flow  and  E^.  grows  linearly,  on  the  average. 

Figure  16  shows  that  breaking  of  cusped  waves 
(upper  dashed  curve  for  Rg  = 200)  is  a more  energetic  process 
than  when  the  waves  do  not  break.  This  is  consistent  with 
an  earlier  observation  that  cusped  waves  are  buoyancy  dominated, 
so  that  they  depend  strongly  on  the  distribution  of  density. 

Figure  17  shows  the  time  dependences  of  the  density 
and  momentum  thicknesses  for  R = 3-5  in  the  cases  J = 0.4  and 
J = 0.1  (Runs  7 and  4).  The  density  thickness  for  the  cusped 
wave  case  (a)  grows  linearly  even  at  very  long  times. 

Figure  18  shows  the  z dependences  of  o,  u and  Rj 
with  R = 3.5>  and  J = 0.4  and  0.1  at  a time  t=130,  compared 
with  their  initial  values  (dashed  lines).  The  strongly 
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unstable  cat's  eye  case  (a,  d and  f)  shows  evidence  for  sheet 
splitting,  characterized  by  a nearly  constant  central  density 
region  bounded  by  regions  of  large  Richardson  number. 

At  the  time  shown  in  Figure  18,  the  average  Richardson 
number  over  the  central  flow  region  (where  R^(z)  r1  0)  is  0.29 
for  the  cat's  eye  case  (f)  and  0.27  for  the  cusped  wave  case 
(e).  This,  together  with  a similar  observation  in  Section  4 , 
is  in  agreement  with  the  experimental  results  of  Thorpe  (1973), 
who  observes  average  gradient  Richardson  numbers  of  order 
0.3  in  the  fully  developed  central  flow  region  resulting 
from  two  dimensional  shear  induced  instabilities.  These  results 
do  not  support  heuristic  arguments  by  Businger  (1969)  that 
the  average  gradient  Richardson  number  should  approach  unity 
at  late  times. 

t 

6 . Conditions  for  Cusped  Wave  Instability  in  the  Malta 

Thermocline 

— 

The  significance  of  cusped  wave  instabilities  in 

I 

the  ocean  is  not  understood.  Since  this  instability  predicts  ;■ 

a slow  but  persistent  linear  spreading  of  the  interfacial 
region,  it  is  of  interest  to  explore  conditions  for  the 
appearance  of  cusped  waves  in  the  thermocline.  In  this 
section  we  examine  how  long  wavelength  internal  waves  are 
modified  by  the  presence  of  a thin  sheet  of  enhanced 
thermal  gradient.  We  will  be  interested  in  conditions  where 
an  appropriately  defined  local  shear  thickness  exceeds  the 
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sheet  thickness  by,  say,  /?.  It  is  useful  to  employ  dimen- 
sional quantities  in  this  discussion. 

We  model  the  mean  thermocline  thermal  gradient 

2 

with  a Gaussian  profile,  as  shown  in  Figure  19.  Here  N , the 

dimensional  mean  Brunt-Vaisala  frequency  in  the  absence  of 

sheets,  is  plotted  versus  depth  in  meters.  The  mean  Malta 

profile  (Woods,  1968)  is  shown  for  comparison.  We  assume 

that  a single  sheet  is  located  at  the  thermocline  center, 

2 

taken  as  z = 0 In  the  following,  and  that  the  total  N is 
given  by 


JT(z)  = 


Ngexp 


- 1 (i)’ 


+ N exp 
m ^ 


- w 


(6.1) 


2 2 

where  N >>  Nn  and  h <<  H.  Values  typical  of  those  in  the 
m u 

2 — 

Malta  thermocline  are  H = 550  cm,  h = 5 cm,  Nq  = 0.0008  sec 
2 _2 

and  Nm  = 0.004  sec  , corresponding  to  total  temperature 
differences  of  6.3°C  and  0.26°C  across  the  thermocline  and 
sheet,  respectively. 

There  are  an  infinite  number  of  internal  wave  modes 

of  given  wavelength  L having  different  real  wave  velocities 

c.  We  have  extensively  studied  properties  of  the  lowest 

shear  mode  (whose  maximum  shear  occurs  at  z = 0)  as  a 

~2 

function  of  L/h,  H/h  and  NQ.  This  involves  solution  of  the 
eigenvalue  equation  for  internal  waves  in  the  Boussinesq 
approximation  (Phillips,  1969,  pg.  1 6 1 ) with  boundaries  at 
z = i 00  (no  surface  waver.)  . The  results  are  summarised  in 
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the  wave  stability  diagrams  of  FI  pi  ire. a 20  and  21.  Hero, 
we  eliminate  the  unknown  velocity  Un  in  the  overall 
Richardson  number  J in  favor  of  the  internal  wave  amplitude 
a and  define  a modified  Richardson  number 

5 - (!)"->•  <6-2> 

We  also  define,  somewhat  arbitrarily , a local  shear  thickness 

ratio  R,  equal  to  the  value  of  /? z/h  where  the  computed 

horizontal  wave  velocity  rises  to  0.5 2 of  its  maximum  value 

(this  reduces  to  the  previous  meaning  for  R in  Section  2 

for  horizontal  wave  velocities  and  sheet  density  of  the 

~2 

form  (2.1)  and  (2.2),  where  6 - /tt/2  Rh)  . For  small  Nq, 

corresponding  to  large  thermal  steps,  Figure  20  shows 'that  the 

internal  waves  are  substantially  altered  near  the  sheet  and 

R is  of  order  unity  or  less  except  for  very  long  wavelengths 

-2 

(hatched  region  to  the  right).  For  somewhat  larger  NQ 
(smaller  thermal  steps)  the  sheet  is  less  effective  in 
altering  the  local  shear  thickness,  and  interesting  parameter 
regions  occur  where  R is  substantially  greater  than  unity 
(Figure  21),  where  we  expect  cusped  wave  instability. 

The  J contours  in  Figures  20  and  21  indicate 
that  shear  instability  leading  to  cat's  eyes  (J  < 1/4)  occurs 
only  for  relatively  large  amplitude  internal  waves.  For 
example  with  = 0.2  (Figure  20),  H/h  = 110  and  L/h  = 150, 
it  Is  necessary  that  a/I.  = O.lfi  to  reduce  .1  to  1/4.  Amplitude- 
wavelength  rat  !<).;  i.f  i lil.;  >>rder  .are  .approach  I up  valie  .;  where 


nonlinear  breaking  or out 
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•s.  Till  a point,  in  further  illustrated 
in  Figure  21  for  the  case  of  a weak  thermal  step.  The  heavy 
solid  curve  represents  the  locus  of  points  where  J = 1/14  and, 
simultaneously,  the  mean-square  wave  slope  is  unity.  Above 
this  curve,  nonlinear  breaking  is  more  likely  than  the  cat's 
eye  shear  instability. 

The  results  in  this  section  are  heuristic  in  the 
sense  that  we  have  not  studied  the  detailed  stability  charac- 
teristics of  internal  wave  profiles  modified  by  the  presence 
of  sheets.  This  would  be  a considerably  more  ambitious 
undertaking  which  is  complicated  by  the  periodic  time  vari- 
ation of  the  actual  internal  wave  shear.  The  results  suggest, 
however,  that  cusped  wave  instability  may  be  a commonly 
occurring  mechanism  which  gives  sheet  spreading  (i.e.,  vertical 
heat  transport)  by  relatively  small  amplitude  internal  waves 
with  J > 1/4. 


7 • Conclusion 

We  have  established  conditions  which  lead  to  non- 
linear cusped  waves  in  layered,  shear  flows.  For  the  model 
error  function  mean  velocity  and  density  profiles  used,  we  find 
that  cusped  wave  instability  occurs  when  the  length  scale  of 
the  shear  layer  exceeds  that  of  the  density  transition  layer 
by  /2  and,  simultaneously,  the  centerline  Richardson  number 
exceeds  1/4.  In  other  cases  instability  develops  into  a 
familiar  cat's  eye  pattern.  The  cusped  wave  Instability  is 
of  special  Interest  in  that  It  continues  to  extract  energy 
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from  thr  moan  shear,  whereas  the  cat 'a  ry  t pattern  is  bounded 
in  its  energy  content  and  vertical  extent.  We  have  presented 
indirect  evidence  that  conditions  for  the  cusped  wave  insta- 
bility may  be  relatively  common  for  realistic  thertnoclines , 
sheet  thicknesses  and  internal  wave  amplitudes. 
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FIGURE  CAPTIONS 


Figure  1.  Stability  Boundary  for  R = 1 and  Neutral  Stability 
Curves  for  Various  R (Error  Function  Profiles) . 

Figure  2.  Stability  Boundaries  for  R = 3.5  (Error  Function 
Profiles) . 

Figure  3.  Vertical  z Dependencies  of  the  Real  and  Imaginary 
Parts  of  the  Most  Rapidly  Growing  Eigen  Functions 
for  Various  R and  J. 

Figure  4.  Constant  Density  Contours  to  t = 80.  R = 1, 

J = 0.1,  a = 0.43,  Re  = 100,  Pr  = 10  (Run  1). 

Figure  5-  Constant  Vorticity  Contours  (Solid)  and 
Streamlines  (Dotted)  Corresponding  to  Figure  4. 

Figure  6.  Vertical  Profiles  of  the  Mean  Density  Perturbation 
o = 1000* (p-l;.  Mean  Velocity  u and  Mean  Richardson  Number 
R’  at  Times  t = 0 (Dashed)  and  150.  R = 1,  J = 0.1, 

a = 0.43,  Re  = 100,  Pr  = 10  (Run  2). 

Figure  7.  Momentum  and  Density  Thicknesses  Au  and  Ap 

(Normalized  by  their  Initial  Values)  Versus  Time. 

Conditions  of  Figure  6. 

Figure  8.  Fluctuation  Energy,  Reynolds  Shear  Stress, 

Buoyancy  Flux  and  Dissipation  Rate  Versus  Time.  Conditions 
Figure  6. 

Figure  9-  Density  Contours  to  t = 80.  R = 3*5,  J = 0.1, 
a = 0.38,  Re  = 100,  Pp  = 10  (Run  3). 

Figure  10.  Vorticity  Contours  (Solid)  and  Streamlines 
(Dotted)  Corresponding  to  Figure  9. 

Figure  11.  Fluctuation  Energy,  Reynolds  Shear  Stress, 

Buoyancy  Flux  and  Dissipation  Rate  Versus  Time.  R = 3-5, 

J = 0.1,  a = 0.38,  Re  = 100,  Pr  = 10  (Run  4). 

Figure  12.  Cusped  Wave  Constant  Density  Contours  for 
R = 3.5,  J = 0.4,  ct  = 0.5,  Re  = 100,  Pr  = 10  (Run  7). 

Figure  13.  Vorticity  Contours  and  Streamlines  Corresponding 
to  Figure  12. 
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Figure  14.  Long  Time  Behavior  Corresponding  to  Figures  IB 
and  13.  The  Arrows  Indicate  the  Position  of  a Secondary 
Wave  Moving  to  the  Left. 

p 

Figure  15.  Breaking  of  Cussed  Wave  Density  Contours  for 
R = 3-5,  J = 0.4,  a = 0.5,  Re  = 200,  Pr  = 10  (Run  9). 

Figure  16.  Fluctuation  Energy,  Reynolds  Shear  Stress, 
Buoyancy  Flux  and  Dissipation  Versus  Time.  R = 3-5, 

J = 0.4,  a = 0.5,  Pr  = 10,  Rg  = 100  (Run  7).  The  Upper 

Dashed  Curve  Shows  at  Rg  = 200  (Run  9)  for  a Breaking 

Cusped  Wave. 
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Figure  17.  Density  and  Momentum  Thicknesses  Versus  Time 
for  R = 3.5  and  ?r  = 10.  (a)  and  (c)  Give  Cusped  Wave 

Formation  with  J = 0.4  (Run  7).  (b)  and  (d)  Give  Cat's 

Eye  Formation  with  J = 0.1  (Run  4). 

Figure  18.  Vertical  Profiles  of  the  Mean  Density  Perturbation 
and  Mean  Velocity  at  Times  t = 0 (Dashed)  and  130. 

Conditions  of  Figure  17. 

Figure  19.  Mean  Thermocline  Model.  The  Square  Brunt-Vaisala 
Frequency  (sec  ) is  shown  Versus  Depth  (meters). 

Figure  20.  Lowest  Shearing  Mode  Internal  Wave  Stability 

Diagram  for  Nq  = 0.2.  Constant  R ( ) and  J( — ) Contours. 

The  Lowest  Shearing  Mode  Does  Not  Exist  in  the  Left  Hatched 
Region.  The  Circle  Corresponds  Approximately  to  Conditions 
Reported  by  Woods  and  Wiley  (1972)  for  Malta. 

-N.  p 

Figure  21.  Internal  Wave  Stability  Diagram  for  Nn  = 0.4. 

Constant  R ( ) and  J ( — ) Contours.  The  heavy  solid 

line  Is  discussed  in  the  text. 
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